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ABSTRACT 



Context. In Cepheids close to the red edge of the classical instability strip, a coupling occurs between the acoustic oscillations and the 
convective motions close to the surface. The best topical models that account for this coupling rely on 1-D time-dependent convection 
(TDC) formulations. However, their intrinsic weakness comes from the large number of unconstrained free parameters entering in the 
description of turbulent convection. 

Aims. We compare two widely used TDC models with the first two-dimensional nonlinear direct numerical simulations (DNS) of the 
convection-pulsation coupling in which the acoustic oscillations are self- sustained by the /c-mechanism. 

Methods. The free parameters appearing in the Stellingwerf and KuhfuB TDC recipes are constrained using a -test with the time- 
dependent convective flux that evolves in nonlinear simulations of highly-compressible convection with /c-mechanism. 
Results. This work emphasises some inherent limits of TDC models, that is, the temporal variability and non-universality of their free 
parameters. More importantly, within these limits, Stellingwerf 's formalism is found to give better spatial and temporal agreements 
with the nonlinear simulation than KuhfuB's one. It may therefore be preferred in 1-D TDC hydrocodes or stellar evolution codes. 
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1. Introduction 

The first theoretical calculations of the Cepheids instability strip 
done in the 60 's assumed that convection was steady with re- 
spect to oscillations. Unfortunately, this "frozen-in convection" 
approximation led to a cooler red edge than the observed one 
as the strong coupling between convection and pulsations oc- 
curring in cool Cepheids was ignored (e.g. Baker & Kippenhahn 
1965). Then, following the pioneering works of Unno (1967) 
and Gough (1977), several time-dependent convection (TDC) 
models have been developed to investigate the influence of 
the convection onto the pulsational stability (e.g. Stellingwerf 
1982; KuhfuB 1986; Gehmeyr & Winkler 1992a). The last up- 
to-date TDC models actually succeed in predicting both a red 
edge close to the observed one and realistic luminosity curves 
(e.g. Yeckoetal. 1998; Bono et al. 1999; Feuchtinger 1999; 
Kollath et al. 2002). 

However, all of these TDC models suflTer from a common 
weakness due to the numerous free parameters, usually known 
as a coeflftcients, that describe the turbulent convection (e.g. the 
8 dimensionless parameters in Smolec & Moskalik 2010). These 
parameters are either obtained from a fit to the observations or 
hardly constrained by theory when taking the asymptotic limit 
of stationary convection (Vitense 1953; Bohm-Vitense 1958). A 
parametric study carried out by Yecko et al. (1998) has empha- 
sised the intrinsic degeneracy of TDC models as similar instabil- 
ity strips have been obtained with diflTerent sets of parameters. 

Direct numerical simulations (hereafter DNS) are able to 
constrain these TDC models as they fully account for the non- 
linearities involved in the convection-pulsation coupling. These 
nonlinear simulations are challenging as they require both large- 
amplitude oscillations and convective motions. In our last 2-D 



simulations, we have improved the way acoustic waves are gen- 
erated by reproducing the self-consistent excitation operating 
in variable stars, that is, the A:-mechanism (Gastine & Dintrans 
2011, hereafter GD2011). The resulting coupling of acoustic 
modes with convection is therefore more consistent as the mode 
amplitude is not imposed artificially. We have shown in GD201 1 
that the convective plumes may either quench the radial oscilla- 
tions or coexist with the acoustic modes, depending mainly on 
the density contrast of the equilibrium model. 

The purpose of this letter is to compare the fully nonlinear 
results with two main TDC models: (/) the first one refers to 
an initial formulation of Stellingwerf (1982) that has been used 
and improved by Bono & Stellingwerf (1994) and Bono et al. 
(1999); (//) the other one has been developed by KuhfuB 
(1986) and Gehmeyr & Winkler (1992a) and is implemented 
both in the Vienna hydrocode (e.g. Wuchterl & Feuchtinger 
1998; Feuchtinger 1999) and in the Florida-Budapest one (e.g 
Yeckoetal. 1998; Kollath etal. 2002). In these models, a sin- 
gle equation for the turbulent kinetic energy St is added to the 
classical mean-field equations and the main second-order corre- 
lations, as for example the convective flux, are expressed as a 
function of St only (see e.g. Baker 1987; Gehmeyr & Winkler 
1992b; Buchler 2000, 2009). 

We investigate here in more details a particular simulation 
where the oscillations strongly modulate the convective flux. The 
nonlinear results are first compared at each snapshot with the 
TDC recipes by computing a ;^^-statistics of the relevant a co- 
eflftcients. Secondly, the mean values of a are used to compare 
the optimal TDC fluxes with the DNS ones. The formulation of 
Stellingwerf is found to be closer to the nonlinear result than 
KuhfuB's one: (/) the temporal variation of the a coefficient is 
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weaker; (//) the mean convective flux is closer to the DNS one, 
especially in the description of the overshooting area. 

2. The hydrodynamic model 

We consider a local 2-D box of size Lx x L^, filled by a per- 
fect monatomic gas, and centered on both sides of an ionisa- 
tion region, of which the associated opacity bump is shaped by a 
hollow in the temperature-dependent radiative conductivity pro- 
file K(T). Such a configuration can lead to unstable acoustic 
modes due to the driving term V • K(T)VT in the energy equa- 
tion, i.e. this is the A:-mechanism (Gastine & Dintrans 2008a,b). 
Furthermore, this hollow in K(T) is deep enough such as the 
equilibrium temperature gradient is locally superadiabatic and 
convective motions develop here according to Schwarzschild's 
criterion. 

The governing hydrodynamic equations are written in nondi- 
mensional form by choosing as the length scale and ^Jc^T\^ 

as the velocity one, hence the time scale L^/ yJc^T^ (with Cp 
the specific heat and Ttop the surface temperature). The result- 
ing nonlinear set of equations is advanced in time with the high- 
order finite-diff'erence pencil code\ which is fully explicit except 
for the radiative diff'usion term that is solved implicitly thanks to 
a parallel alternate direction implicit (ADI) solver. With the cho- 
sen units, the simulation box spans about 10% of the star radius 
on both sides of the ionisation region while the timestep is about 
1 minute, such that a simulation typically spans over 4500 days 
(see GD2011). 

3. Results 

The time evolution of the convective flux obtained in fully non- 
linear 2-D simulations is compared with the following TDC ex- 
pressions developed by Stellingwerf (1982) and KuhfuB (1986): 



rst(z, = ast-St sign(V - Vad) VlV - Va, 

^ TKuiz, t) = aKuA ^|St(V- Vad) , 
where V = Jin T/dln p, Vad = 1 - (^v/(^p and 



(1) 



A = cp <p) (T) and B = (T) Vad, (2) 



with p and p the pressure and density, respectively, and the 
brackets denote an horizontal average. Two free dimensionless 
parameters Ofst and a^u are introduced that the simulations al- 
low to constrain. We first focus on a 2-D DNS that is similar 
to the G8 simulation in GD2011, that is, a simulation in which 
the total kinetic energy is almost entirely contained in acoustic 
modes excited by /c-mechanism (80%), the rest being in convec- 
tive plumes (20%). The convective motions do not aff'ect much 
the growth and the nonlinear saturation of the unstable radial 
acoustic mode. We can therefore expect that the heat transport is 
strongly modulated by wave motions, what is an ideal frame for 
testing the accuracy of time-dependent convection models. 

3.1. Temporal modulation of the convective flux in the DNS 

The imposed bottom flux Fbot is mainly transported through the 
computational domain by the radiative Trad, enthalpy and kinetic 

^ See http://www.nordita.org/software/pencil-code/ and 
Brandenburg & Dobler (2002). 
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Fig. 1. a) Mean vertical profiles of radiative Trad (solid blue line), 
turbulent enthalpy Tcom (dashed green line), kinetic Tkin (dotted 
black line), modes Tmod (dot-dashed magenta line) and total Ttot 
(dotted red line) fluxes, normalised to the bottom flux Fbot- b) 
Temporal power spectrum for the convective flux only. 



Tkin fluxes. Following GD201 1, the enthalpy flux is divided into 
the classical convective flux Tcom and the contribution coming 
from the (unstable) fundamental mode (hereafter Tmod) with 



Tconyiz, t) = Cp {puj') , ?mod(^, = Cp {pu^} 0, 

rradiz, t) =- {K(T)VT} , T^niz, t)=^ {pu,u^) , 



(3) 



where the primed quantities denote the fluctuations about the 
horizontal average and 6 is the temperature eigenfunction of the 
fundamental mode. The resulting time-averaged and normalised 
fluxes are given in Fig. la. 

The bulk of the total flux is transported by the radiative flux, 
except in the convective zone where ?^onv/^bot - 20%, while the 
kinetic flux is negligible (^Fkin/^bot ^ !%)• Concerning Tmod, 
one notes that it is hardly as large as Tcorw in the convective 
zone. This quantity is a good signature of the amplitude of the 
acoustic modes as the higher Tmod, the larger the radial oscil- 
lations (GD2011), therefore confirming the efficiency of the k- 
mechanism in this simulation. 

The signature of the temporal modulation of the convective 
flux is extracted from its Fourier spectrum in time, that is, we 
first compute Tcorwiz, oj), with oj the frequency, and second in- 
tegrate over the vertical direction z to get the power spectrum 
Tcorwi^) (Fig. lb). Several discrete peaks appear about given 
frequencies, of which the physical nature is emphasised after 
superimposing the linear acoustic eigenfrequencies (the vertical 
dashed blue lines). The matching with the fundamental mode 
frequency ojqo = 3.85 is perfect while the weak- amplitude 
secondary peaks rather correspond to harmonics of ojqo (i.e. 
2ajoo, 3ojoo, . . . , shaped by downward-directed vertical gray ar- 
rows in Fig. lb) than to the acoustic overtones ojqi, ojo2, It 

means that the amplitude of the fundamental acoustic mode is 
large enough to generate several harmonics through a nonlinear 
cascade and this is again an indication of both the robustness 
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Fig. 2. Evolution of the convective flux in a (t, z) plane: (a) in the DNS; (b) with Stellingwerf 's formalism 9^su (c) with KuhfuB's for- 
malism Tku (for ast = o;ku = !)• The vertical extent is centered on both sides of the convection zone to emphasise its oscillations. 
These snapshots are computed after 1800 periods of oscillations, i.e. well after the nonlinear saturation is achieved. 



of the /c-driving and the relevance of this kind of convection- 
pulsation simulation to check the TDC recipes. 

3.2. DNS vs TDC models: convective patterns 

We first compare the time evolution of the horizontally aver- 
aged convective flux obtained in the DNS (Eq. 3) with its TDC 
counterparts Tst and Tku (Eq. 1). As we are interested here in 
the qualitative agreement between the convective patterns in the 
plane (t, z), we simply assume a^t = c^Ku = 1- 

The three resulting patterns are displayed in Fig. 2 for a time 
interval spanning about 4 periods of the fundamental acoustic 
mode. The black areas denote positive values for the convective 
flux and therefore delimit the convective zone. An oscillatory be- 
haviour is clearly visible in each panel, with a period that looks 
similar to the one of the fundamental acoustic mode, that is, al- 
most 4 oscillation cycles are depicted. This is consistent with the 
large peak shown around ojoq in the power spectrum in Fig. lb. 

The two TDC patterns also display a good agreement 
with the nonlinear simulation regarding the overshooting phe- 
nomenon. We indeed recover the same dark red structures below 
the convective zone that correspond to the downdrafts entering 
in the radiative zone (this penetration being associated with a 
negative convective flux). Nevertheless, we note that the over- 
shooting seems to be more vigorous in KuhfuB's model than in 
Stellingwerf 's one as these dark red structures fill in Fig. 2c a 
larger surface in the bottom radiative zone than in Fig. 2b. 

3.3. DNS vs TDC models: statistics of coefficients a 

A one-to-one comparison between the convective flux in the 
DNS and its TDC predictions requires to find the optimal val- 
ues of ast and ^ku- This is done by performing several -tests 
at diff'erent snapshots in the simulation to track the variations of 
coefficients a. The resulting fluctuations of a^t and q^ku around 
their mean values are shown in Fig. 3. 

The dispersion of the Stellingwerf coefficient a^st (upper 
panel) is weaker than the KuhfuB one (lower panel) as its val- 
ues are almost within a 5% range around the mean osi = 0.462. 
On the contrary, several outliers with values greater than 1 0% are 
found in the evolution of a^u which then appears more chaotic 
around the mean value a^u = 0.076. As a consequence, the rel- 
ative standard deviation (depicted in gray in Fig. 3) is weaker in 
Stellingwerf 's case than in KuhfuB's one. 
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Fig. 3. Time evolution around the mean value of coefficients ast 
(upper panel) and o^ku (bottom panel). Horizontal gray spans 
mark the limits of the associated relative standard deviation. 



3.4. DNS vs TDC models: convective fluxes with optimal a's 

We recall that TDC models assume that the coeflftcients a enter- 
ing in Eq. 1 are constant. By adjusting the TDC recipes with the 
instantaneous convective flux throughout the nonlinear simula- 
tion, the optimal value of these coefficients has been deduced. 
The final test, given in Fig. 4, then consists in the comparison 
between the mean nonlinear convective flux taken over the en- 
tire simulation and its best TDC approximations built from these 
optimal a values. 

This figure emphasises that Stellingwerf 's formulation gives 
a better agreement than KuhfuB's one. Indeed, the KuhfuB model 
overestimates the overshooting as the (negative) convective flux 
remains non-negligible until the bottom of the radiative zone. 
On the contrary, the Stellingwerf profile better accounts for 
the local penetration of convective plumes and shows the same 
exponential-like decay of the negative convective flux when 
sinking in the radiative zone (e.g. Dintrans 2009). However, the 
two models are rather similar in the bulk of the convective zone 
where convection is fully developed. One also notes that they 
both predict a negative flux at the top of the convective zone, 
that is, an upper overshooting of convection motions near the 
surface that is not observed in the nonlinear simulation. 
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Fig. 4. Mean convective flux in the DNS (solid black line), com- 
pared with the best TDC predictions based on the optimal a val- 
ues that came out of the -test in §3.3 (Stellingwerf 's model in 
dashed blue line and KuhfuB's one in dotted green line). 



4. Conclusion 

The main weakness of theories of time-dependent convection 
(TDC) lies in the large number of free parameters. This is partic- 
ularly awkward when convection is strongly coupled with pulsa- 
tions as, for example, near the red border of the Cepheid insta- 
bility strip where similar results are obtained with difl'erent sets 
of parameters (Yecko et al. 1998). Additional constraints must 
be found to reduce the intrinsic degeneracy of models. But the 
modelling of the interplay between the turbulent convection and 
oscillations is really a difficult task due to the strong nonlineari- 
ties involved in this coupling. 

One solution to tackle this problem and to bring new con- 
straints consists in performing fully nonlinear simulations, and 
this is the path we have chosen in this work following our first 
study in GD2011. Two widely used TDC theories, namely the 
Stellingwerf (1982) and KuhfuB (1986) ones, are compared with 
results coming from 2-D nonlinear simulations of compress- 
ible convection in which strong acoustic oscillations are self- 
sustained by the A:-mechanism. The heat transport is then mod- 
ulated by the fundamental acoustic mode such that this kind 
of simulation is relevant to investigate the convection-pulsation 
coupling (Figs. 1-2). 

Focusing on the two TDC formulae for the convective flux, 
we compute the evolution of free coefficients "a" from a;^^-test 
applied to the fully nonlinear results (Fig. 3). A large temporal 
variability is found in both cases that weakens the robustness of 
the TDC assumption a = const. Moreover, the mean values a 
are not universal. By applying the same method to other sim- 
ulations performed in GD2011 (Table 1), we indeed do not re- 
cover the same a with a relative standard deviation of about 12% 
(Stellingwerf) and 18% (KuhfuB). 

Within these limits, Stellingwerf 's formulation is found to 
give a better agreement with the nonlinear simulations than 
KuhfuB's one: (/) the final mean convective flux Tst is closer 
to its DNS counterpart, with a much better estimation of the bot- 
tom overshooting; (//) the temporal dispersion of the Ofst coeflft- 
cient is weaker, then its enhanced stability (Fig. 4). This result 
means that the time- dependent convective flux better scales with 
a law ^Tconv ^ St VV - Vad than Tcorw ^ VSI^(V - Vad), and the 



Table 1. Values of the optimal Stellingwerf and KuhfuB a coef- 
ficients pulled out from the nonlinear simulations in GD201 1 . 



Simulation 






G8 


0.46 


0.076 


G8H9 


0.47 


0.099 


G8H8 


0.33 


0.113 


G7 


0.38 


0.098 


G6 


0.38 


0.102 


G6F7 


0.40 


0.082 


G6F5 


0.38 


0.067 



Stellingwerf formulation may probably be preferred in the 1 -D 
hydrocode used in, e.g., the topical Cepheids models. However, 
this study emphasises that both formalisms lead to an artificial 
overshooting at the top of the convection zone. We also note that 
this TDC test involves the exact value of the turbulent kinetic 
energy St provided by the nonlinear simulation, and not by the 
dedicated 1-D TDC equation which is inherently less accurate. 
As a consequence, the obtained profiles are certainly the best we 
can expect from these TDC recipes. 

In this work, the temporal modulation of the convective flux 
is ensured by the acoustic modes excited by A:-mechanism. An 
interesting prospect could be the case of a modulation based on 
the internal gravity waves excited by convection itself. Indeed, 
convection can excite gravity waves in variable stars, either by 
the means of the penetration of convective elements into stably 
stratified regions as in solar-type stars (e.g. Dintrans et al. 2005), 
or through the so-called "convective blocking" mechanism as in 
white dwarfs or Gamma Doradus stars (e.g. Pesnell 1987). In 
both cases, the convective flux is ipso facto modulated by gravity 
waves and it may be interesting in that respect to also check the 
accuracy of time-dependent convection models. 
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Appendix A: Setup of the simulations 

A.1. The equilibrium model 

As in GD201 1, our system represents a zoom around an ionisa- 
tion region. As we are computing local simulations, the vertical 
gravity g = -ge^ and the kinematic viscosity v are assumed 
to be constant. Following our purely radiative model of the k- 
mechanism (Gastine & Dintrans 2008a,b), the ionisation region 
is represented by a temperature-dependent radiative conductivity 
profile that mimics an opacity bump: 



K(T) = 



with 



Km 



1 +^ 



-7r/2-harctan(crr+r-) 



7t/2-\- arctan((T^2) 



T- = T-T] 



bump 



(A.1) 



(A.2) 



where Tbump is the position of the hollow in temperature and cr, 
e, and ^ denote its slope, width, and relative amplitude, respec- 
tively. We assume both radiative and hydrostatic equilibria; that 



IS, 



dpo 
dz 

dTo 
dz 



(A3) 



where Fbot is the imposed bottom flux. Following GD2011, we 
chose Lz as the length scale, i.e. [x] = L^, top density ptop and top 
temperature Ttop as density and temperature scales, respectively. 
The velocity scale is then ^Jc^J\^, while time is given in units 

of [t] = LJ ^JcpTtop. 

Table A. 1 then summarises in these dimensionless units the 
parameters of the numerical simulations presented in this study. 
The penultimate column of this table contains the value of the 
frequency ojqo of the fundamental unstable radial mode excited 
by the /c-mechanism, which lies between 3 and 4 for every DNS. 
The last column gives the value of the Rayleigh number, which 
quantifies the strength of the convective motions. It is given by 



where p, m, and T denote density, velocity, and temperature, 
respectively, while K(T) is given by Eq. A.l. The operator 
D/Dt = d/dt -\- u ' V is the usual total derivative, while S is 
the (traceless) rate-of- strain tensor given by 



1 / dui duj 



dxi 



2 \dxj 



divM 

3 ^ 



(A.6) 



Finally, we impose that all fields are periodic in the horizontal 
direction, while stress-free boundary conditions (i.e. = and 
dux/dz = 0) are assumed for the velocity in the vertical one. 
Concerning the temperature, a perfect conductor at the bottom 
(i.e. flux imposed) and a perfect insulator at the top (i.e. temper- 
ature imposed) are applied. 

In order to ensure that both the nonlinear saturation and ther- 
mal relaxation are achieved, the simulations were computed over 
very long times, typically t > 3000. As the eigenfrequency of the 
unstable acoustic mode cooo e [3-4] (see Table A.l), this corre- 
sponds approximately to 1800 periods of oscillation. 



Ra: 



yxcp 



(A.4) 



where Lconv is the width of the convective zone,;^ = Ko/poCp the 
radiative difl'usivity, and s the entropy. 



A.2. The nonlinear equations 

With the parallel alternate direction implicit solver for the radia- 
tive diffusion implemented in the pencil code (see GD201 1), we 
advance the following hydrodynamic equations in time: 

( Dlnp 

= - dlVM, 

Dt 

Du 1 

— - -V/?-h^-h2y(V-S + S- Vlnp), (A.5) 



= — 6iNK{T)VT - (r - l)rdivM + 2pyS^ 

Dt pCy 
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Table A.l. Dimensionless parameters of the numerical simulations. 



DNS 


Gravity g 


Flux Fbot 


Conductivity /^max 




Width e 


Slope cr 


Viscosity v 


Frequency 6l>oo 


Rayleigh 


G8 


8 


4.5 X 10-2 


10-2 


6 


1 


1 


2.5 X 10-4 


3.85 


105 


G8H9 


8 


4.5 X 10-2 


9 X 10-3 


7 


1 


1 


5 X 10-4 


3.83 




G8H8 


8 


4.5 X 10-2 


8 X 10-3 


7.5 


1 


1.1 


5 X 10-4 


3.80 


2x10' 


G7 


7 


4 X 10-2 


10-2 


5.5 


1.5 


0.8 


2.5 X 10-4 


3.62 


8x 104 


G6 


6 


4 X 10-2 


10-2 


6 


1.5 


0.8 


5 X 10-4 


3.35 


8x 104 


G6F7 


6 


3.7 X 10-2 


10-2 


5.7 


1.5 


0.8 


3.5 X 10-4 


3.36 


9x 104 


G6F5 


6 


3.5 X 10-2 


10-2 


5.5 


1.5 


0.8 


3 X 10-4 


3.36 


9x 104 



Notes. The bold-typed one emphasises the DNS mainly discussed in this study. For all these simulations, we assume Ttop = 2, ptop = 10 2, 
= 0.7, and LJL, = 4. 



